/* -*- mode: C -*-  */
/* 
   IGraph library.
   Copyright (C) 2003-2012  Gabor Csardi <csardi.gabor@gmail.com>
   334 Harvard street, Cambridge, MA 02139 USA
   
   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2 of the License, or
   (at your option) any later version.
   
   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.
   
   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software
   Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 
   02110-1301 USA

*/

#include "igraph_memory.h"
#include "igraph_random.h"
#include "igraph_error.h"

#include <assert.h>
#include <string.h> 		/* memcpy & co. */
#include <stdlib.h>

/** 
 * \section about_igraph_matrix_t_objects About \type igraph_matrix_t objects
 * 
 * <para>This type is just an interface to \type igraph_vector_t.</para>
 *
 * <para>The \type igraph_matrix_t type usually stores n
 * elements in O(n) space, but not always. See the documentation of
 * the vector type.</para>
 */

/**
 * \section igraph_matrix_constructor_and_destructor Matrix constructors and
 * destructors
 */

/**
 * \ingroup matrix
 * \function igraph_matrix_init
 * \brief Initializes a matrix.
 * 
 * </para><para>
 * Every matrix needs to be initialized before using it. This is done
 * by calling this function. A matrix has to be destroyed if it is not
 * needed any more; see \ref igraph_matrix_destroy().
 * \param m Pointer to a not yet initialized matrix object to be
 *        initialized. 
 * \param nrow The number of rows in the matrix.
 * \param ncol The number of columns in the matrix.
 * \return Error code.
 *
 * Time complexity: usually O(n), 
 * n is the
 * number of elements in the matrix.
 */

int FUNCTION(igraph_matrix,init)(TYPE(igraph_matrix) *m, long int nrow, long int ncol) {
  int ret1;
  ret1=FUNCTION(igraph_vector,init)(&m->data, nrow*ncol);
  m->nrow=nrow;
  m->ncol=ncol;
  return ret1;
}

/** 
 * \ingroup matrix
 * \function igraph_matrix_destroy
 * \brief Destroys a matrix object.
 * 
 * </para><para>
 * This function frees all the memory allocated for a matrix
 * object. The destroyed object needs to be reinitialized before using
 * it again.
 * \param m The matrix to destroy.
 * 
 * Time complexity: operating system dependent.
 */ 

void FUNCTION(igraph_matrix,destroy)(TYPE(igraph_matrix) *m) {
  FUNCTION(igraph_vector,destroy)(&m->data);
}

/**
 * \ingroup matrix
 * \function igraph_matrix_capacity
 * \brief Returns the number of elements allocated for a matrix.
 * 
 * Note that this might be different from the size of the matrix (as 
 * queried by \ref igraph_matrix_size(), and specifies how many elements 
 * the matrix can hold, without reallocation.
 * \param v Pointer to the (previously initialized) matrix object
 *          to query.
 * \return The allocated capacity.
 * 
 * \sa \ref igraph_matrix_size(), \ref igraph_matrix_nrow(),
 * \ref igraph_matrix_ncol(). 
 * 
 * Time complexity: O(1).
 */

long int FUNCTION(igraph_matrix,capacity)(const TYPE(igraph_matrix) *m) {
  return FUNCTION(igraph_vector,capacity)(&m->data);
}


/**
 * \section igraph_matrix_accessing_elements Accessing elements of a matrix
 */

/**
 * \ingroup matrix
 * \function igraph_matrix_resize
 * \brief Resizes a matrix.
 *
 * </para><para>
 * This function resizes a matrix by adding more elements to it.
 * The matrix contains arbitrary data after resizing it.
 * That is, after calling this function you cannot expect that element
 * (i,j) in the matrix remains the
 * same as before.  
 * \param m Pointer to an already initialized matrix object.
 * \param nrow The number of rows in the resized matrix.
 * \param ncol The number of columns in the resized matrix.
 * \return Error code.
 * 
 * Time complexity: O(1) if the
 * matrix gets smaller, usually O(n)
 * if it gets larger, n is the 
 * number of elements in the resized matrix.
 */

int FUNCTION(igraph_matrix,resize)(TYPE(igraph_matrix) *m, long int nrow, long int ncol) {
  FUNCTION(igraph_vector,resize)(&m->data, nrow*ncol);
  m->nrow=nrow;
  m->ncol=ncol;
  return 0;
}

/**
 * \ingroup matrix
 * \function igraph_matrix_resize_min
 * \brief Deallocates unused memory for a matrix.
 * 
 * </para><para>
 * Note that this function might fail if there is not enough memory
 * available. 
 * 
 * </para><para>
 * Also note, that this function leaves the matrix intact, i.e. 
 * it does not destroy any of the elements. However, usually it involves 
 * copying the matrix in memory.
 * \param m Pointer to an initialized matrix.
 * \return Error code.
 * 
 * \sa \ref igraph_matrix_resize().
 * 
 * Time complexity: operating system dependent.
 */

int FUNCTION(igraph_matrix,resize_min)(TYPE(igraph_matrix) *m) {
  TYPE(igraph_vector) tmp;
  long int size=FUNCTION(igraph_matrix,size)(m);
  long int capacity=FUNCTION(igraph_matrix,capacity)(m);
  if (size == capacity) { return 0; }
  
  IGRAPH_CHECK(FUNCTION(igraph_vector,init)(&tmp, size));
  FUNCTION(igraph_vector,update)(&tmp, &m->data);
  FUNCTION(igraph_vector,destroy)(&m->data);
  m->data = tmp;
  
  return 0;
}


/**
 * \ingroup matrix
 * \function igraph_matrix_size
 * \brief The number of elements in a matrix.
 * 
 * \param m Pointer to an initialized matrix object.
 * \return The size of the matrix.
 *
 * Time complexity: O(1).
 */

long int FUNCTION(igraph_matrix,size)(const TYPE(igraph_matrix) *m) {
  return (m->nrow) * (m->ncol);
}

/**
 * \ingroup matrix
 * \function igraph_matrix_nrow
 * \brief The number of rows in a matrix.
 * 
 * \param m Pointer to an initialized matrix object.
 * \return The number of rows in the matrix.
 * 
 * Time complexity: O(1).
 */

long int FUNCTION(igraph_matrix,nrow)(const TYPE(igraph_matrix) *m) {
  return m->nrow;
}

/**
 * \ingroup matrix
 * \function igraph_matrix_ncol
 * \brief The number of columns in a matrix.
 * 
 * \param m Pointer to an initialized matrix object.
 * \return The number of columns in the matrix.
 * 
 * Time complexity: O(1).
 */

long int FUNCTION(igraph_matrix,ncol)(const TYPE(igraph_matrix) *m) {
  return m->ncol;
}

/** 
 * \ingroup matrix
 * \function igraph_matrix_copy_to
 * \brief Copies a matrix to a regular C array.
 *
 * </para><para>
 * The matrix is copied columnwise, as this is the format most
 * programs and languages use.
 * The C array should be of sufficient size; there are (of course) no
 * range checks.
 * \param m Pointer to an initialized matrix object.
 * \param to Pointer to a C array; the place to copy the data to.
 * \return Error code.
 *
 * Time complexity: O(n),
 * n is the number of 
 * elements in the matrix.
 */

void FUNCTION(igraph_matrix,copy_to)(const TYPE(igraph_matrix) *m, BASE *to) {
  FUNCTION(igraph_vector,copy_to)(&m->data, to);
}

/** 
 * \ingroup matrix
 * \function igraph_matrix_null
 * \brief Sets all elements in a matrix to zero.
 * 
 * \param m Pointer to an initialized matrix object.
 * 
 * Time complexity: O(n),
 * n is the number of  elements in
 * the matrix. 
 */

void FUNCTION(igraph_matrix,null)(TYPE(igraph_matrix) *m) {
  FUNCTION(igraph_vector,null)(&m->data);
}

/**
 * \ingroup matrix
 * \function igraph_matrix_add_cols
 * \brief Adds columns to a matrix.
 * \param m The matrix object.
 * \param n The number of columns to add.
 * \return Error code, \c IGRAPH_ENOMEM if there is
 *   not enough memory to perform the operation.
 *
 * Time complexity: linear with the number of elements of the new,
 * resized matrix.
 */

int FUNCTION(igraph_matrix,add_cols)(TYPE(igraph_matrix) *m, long int n) {
  FUNCTION(igraph_matrix,resize)(m, m->nrow, m->ncol+n);
  return 0;
}

/**
 * \ingroup matrix
 * \function igraph_matrix_add_rows
 * \brief Adds rows to a matrix.
 * \param m The matrix object.
 * \param n The number of rows to add.
 * \return Error code, \c IGRAPH_ENOMEM if there
 *   isn't enough memory for the operation.
 * 
 * Time complexity: linear with the number of elements of the new,
 * resized matrix.
 */

int FUNCTION(igraph_matrix,add_rows)(TYPE(igraph_matrix) *m, long int n) {
  long int i;
  FUNCTION(igraph_vector,resize)(&m->data, (m->ncol)*(m->nrow+n));
  for (i=m->ncol-1; i>=0; i--) {
    FUNCTION(igraph_vector,move_interval2)(&m->data, (m->nrow)*i, (m->nrow)*(i+1),
					  (m->nrow+n)*i);
  }
  m->nrow += n;
  return 0;
}

/**
 * \ingroup matrix
 * \function igraph_matrix_remove_col
 * \brief Removes a column from a matrix.
 * 
 * \param m The matrix object.
 * \param col The column to remove.
 * \return Error code, always returns with success. 
 * 
 * Time complexity: linear with the number of elements of the new,
 * resized matrix.
 */

int FUNCTION(igraph_matrix,remove_col)(TYPE(igraph_matrix) *m, long int col) {
  FUNCTION(igraph_vector,remove_section)(&m->data, (m->nrow)*col, (m->nrow)*(col+1));
  m->ncol--;
  return 0;
}

/**
 * \ingroup matrix
 * \function igraph_matrix_permdelete_rows
 * \brief Removes rows from a matrix (for internal use).
 * 
 * Time complexity: linear with the number of elements of the original
 * matrix. 
 */

int FUNCTION(igraph_matrix,permdelete_rows)(TYPE(igraph_matrix) *m, long int *index, long int nremove) {
  long int i, j;
  for (j=0; j<m->nrow; j++) {
    if (index[j] != 0) {
      for (i=0; i<m->ncol; i++) {
        MATRIX(*m, index[j]-1, i) = MATRIX(*m, j, i);
      }
    }
  }
  /* Remove unnecessary elements from the end of each column */
  for (i=0; i<m->ncol; i++)
    FUNCTION(igraph_vector,remove_section)(&m->data,
      (i+1)*(m->nrow-nremove), (i+1)*(m->nrow-nremove)+nremove);
  FUNCTION(igraph_matrix,resize)(m, m->nrow-nremove, m->ncol);

  return 0;
}

/**
 * \ingroup matrix
 * \function igraph_matrix_delete_rows_neg
 * \brief Removes columns from a matrix (for internal use).
 * 
 * Time complexity: linear with the number of elements of the original
 * matrix. 
 */

int FUNCTION(igraph_matrix,delete_rows_neg)(TYPE(igraph_matrix) *m, 
					    const igraph_vector_t *neg, long int nremove) {
  long int i, j, idx=0;
  for (i=0; i<m->ncol; i++) {
    for (j=0; j<m->nrow; j++) {
      if (VECTOR(*neg)[j] >= 0) {
	MATRIX(*m, idx++, i) = MATRIX(*m, j, i);
      } 
    }
    idx=0;
  }
  FUNCTION(igraph_matrix,resize)(m, m->nrow-nremove, m->ncol);

  return 0;
}

/**
 * \ingroup matrix
 * \function igraph_matrix_copy
 * \brief Copies a matrix.
 *
 * </para><para>
 * Creates a matrix object by copying from an existing matrix.
 * \param to Pointer to an uninitialized matrix object.
 * \param from The initialized matrix object to copy.
 * \return Error code, \c IGRAPH_ENOMEM if there
 *   isn't enough memory to allocate the new matrix.
 * 
 * Time complexity: O(n), the number
 * of elements in the matrix.
 */

int FUNCTION(igraph_matrix,copy)(TYPE(igraph_matrix) *to, const TYPE(igraph_matrix) *from) {
  to->nrow = from->nrow;
  to->ncol = from->ncol;
  return FUNCTION(igraph_vector,copy)(&to->data, &from->data);
}

#ifndef NOTORDERED

/**
 * \function igraph_matrix_max
 * 
 * Returns the maximal element of a matrix.
 * \param m The matrix object.
 * \return The maximum element. For empty matrix the returned value is
 * undefined. 
 * 
 * Added in version 0.2.</para><para>
 *
 * Time complexity: O(n), the number of elements in the matrix.
 */

igraph_real_t FUNCTION(igraph_matrix,max)(const TYPE(igraph_matrix) *m) {
  return FUNCTION(igraph_vector,max)(&m->data);
}

#endif

/**
 * \function igraph_matrix_scale
 * 
 * Multiplies each element of the matrix by a constant.
 * \param m The matrix.
 * \param by The constant.
 *
 * Added in version 0.2.</para><para>
 * 
 * Time complexity: O(n), the number of elements in the matrix.
 */

void FUNCTION(igraph_matrix,scale)(TYPE(igraph_matrix) *m, BASE by) {
  FUNCTION(igraph_vector,scale)(&m->data, by);
}

/**
 * \function igraph_matrix_select_rows
 * \brief Select some rows of a matrix.
 * 
 * This function selects some rows of a matrix and returns them in a
 * new matrix. The result matrix should be initialized before calling
 * the function.
 * \param m The input matrix.
 * \param res The result matrix. It should be initialized and will be
 *    resized as needed.
 * \param rows Vector; it contains the row indices (starting with
 *    zero) to extract. Note that no range checking is performed.
 * \return Error code.
 * 
 * Time complexity: O(nm), n is the number of rows, m the number of
 * columns of the result matrix.
 */

int FUNCTION(igraph_matrix,select_rows)(const TYPE(igraph_matrix) *m, 
					TYPE(igraph_matrix) *res, 
					const igraph_vector_t *rows) {
  long int norows=igraph_vector_size(rows);
  long int i, j, ncols=FUNCTION(igraph_matrix,ncol)(m);
  
  IGRAPH_CHECK(FUNCTION(igraph_matrix,resize)(res, norows, ncols));
  for (i=0; i<norows; i++) {
    for (j=0; j<ncols; j++) {
      MATRIX(*res, i, j) = MATRIX(*m, (long int)VECTOR(*rows)[i], j);
    }
  }
  
  return 0;
}

/**
 * \function igraph_matrix_select_rows_cols
 * \brief Select some rows and columns of a matrix.
 * 
 * This function selects some rows and columns of a matrix and returns
 * them in a new matrix. The result matrix should be initialized before
 * calling the function.
 * \param m The input matrix.
 * \param res The result matrix. It should be initialized and will be
 *    resized as needed.
 * \param rows Vector; it contains the row indices (starting with
 *    zero) to extract. Note that no range checking is performed.
 * \param cols Vector; it contains the column indices (starting with
 *    zero) to extract. Note that no range checking is performed.
 * \return Error code.
 * 
 * Time complexity: O(nm), n is the number of rows, m the number of
 * columns of the result matrix.
 */

int FUNCTION(igraph_matrix,select_rows_cols)(const TYPE(igraph_matrix) *m, 
					TYPE(igraph_matrix) *res, 
					const igraph_vector_t *rows,
					const igraph_vector_t *cols) {
  long int nrows=igraph_vector_size(rows);
  long int ncols=igraph_vector_size(cols);
  long int i, j;
  
  IGRAPH_CHECK(FUNCTION(igraph_matrix,resize)(res, nrows, ncols));
  for (i=0; i<nrows; i++) {
    for (j=0; j<ncols; j++) {
      MATRIX(*res, i, j) = MATRIX(*m, (long int)VECTOR(*rows)[i],
	                                  (long int)VECTOR(*cols)[j]);
    }
  }
  
  return 0;
}

/**
 * \function igraph_matrix_get_col
 * \brief Select a column.
 * 
 * Extract a column of a matrix and return it as a vector.
 * \param m The input matrix.
 * \param res The result will we stored in this vector. It should be
 *   initialized and will be resized as needed.
 * \param index The index of the column to select.
 * \return Error code.
 * 
 * Time complexity: O(n), the number of rows in the matrix.
 */

int FUNCTION(igraph_matrix,get_col)(const TYPE(igraph_matrix) *m, 
				    TYPE(igraph_vector) *res,
				    long int index) {
  long int nrow=FUNCTION(igraph_matrix,nrow)(m);

  if (index>=m->ncol) {
    IGRAPH_ERROR("Index out of range for selecting matrix column", IGRAPH_EINVAL);
  }
  IGRAPH_CHECK(FUNCTION(igraph_vector,get_interval)(&m->data, res, 
						    nrow*index, nrow*(index+1)));
  return 0;
}

/**
 * \function igraph_matrix_sum
 * \brief Sum of elements.
 * 
 * Returns the sum of the elements of a matrix.
 * \param m The input matrix.
 * \return The sum of the elements.
 * 
 * Time complexity: O(mn), the number of elements in the matrix.
 */

BASE FUNCTION(igraph_matrix,sum)(const TYPE(igraph_matrix) *m) {
  return FUNCTION(igraph_vector,sum)(&m->data);
}

/**
 * \function igraph_matrix_all_e
 * \brief Are all elements equal?
 * 
 * \param lhs The first matrix.
 * \param rhs The second matrix.
 * \return Positive integer (=true) if the elements in the \p lhs are all 
 *    equal to the corresponding elements in \p rhs. Returns \c 0
 *    (=false) if the dimensions of the matrices don't match.
 * 
 * Time complexity: O(nm), the size of the matrices.
 */

igraph_bool_t FUNCTION(igraph_matrix,all_e)(const TYPE(igraph_matrix) *lhs, 
					    const TYPE(igraph_matrix) *rhs) {
  return lhs->ncol==rhs->ncol && lhs->nrow==rhs->nrow &&
    FUNCTION(igraph_vector,all_e)(&lhs->data, &rhs->data);
}

igraph_bool_t 
FUNCTION(igraph_matrix,is_equal)(const TYPE(igraph_matrix) *lhs, 
				 const TYPE(igraph_matrix) *rhs) {
  return FUNCTION(igraph_matrix,all_e)(lhs, rhs);
}

#ifndef NOTORDERED

/**
 * \function igraph_matrix_all_l
 * \brief Are all elements less?
 * 
 * \param lhs The first matrix.
 * \param rhs The second matrix.
 * \return Positive integer (=true) if the elements in the \p lhs are all 
 *    less than the corresponding elements in \p rhs. Returns \c 0
 *    (=false) if the dimensions of the matrices don't match.
 * 
 * Time complexity: O(nm), the size of the matrices.
 */

igraph_bool_t FUNCTION(igraph_matrix,all_l)(const TYPE(igraph_matrix) *lhs, 
					    const TYPE(igraph_matrix) *rhs) {
  return lhs->ncol==rhs->ncol && lhs->nrow==rhs->nrow &&
    FUNCTION(igraph_vector,all_l)(&lhs->data, &rhs->data);
}

/**
 * \function igraph_matrix_all_g
 * \brief Are all elements greater?
 * 
 * \param lhs The first matrix.
 * \param rhs The second matrix.
 * \return Positive integer (=true) if the elements in the \p lhs are all 
 *    greater than the corresponding elements in \p rhs. Returns \c 0
 *    (=false) if the dimensions of the matrices don't match.
 * 
 * Time complexity: O(nm), the size of the matrices.
 */

igraph_bool_t FUNCTION(igraph_matrix,all_g)(const TYPE(igraph_matrix) *lhs, 
					    const TYPE(igraph_matrix) *rhs) {
  return lhs->ncol==rhs->ncol && lhs->nrow==rhs->nrow &&
    FUNCTION(igraph_vector,all_g)(&lhs->data, &rhs->data);
}

/**
 * \function igraph_matrix_all_le
 * \brief Are all elements less or equal?
 * 
 * \param lhs The first matrix.
 * \param rhs The second matrix.
 * \return Positive integer (=true) if the elements in the \p lhs are all 
 *    less than or equal to the corresponding elements in \p
 *    rhs. Returns \c 0 (=false) if the dimensions of the matrices
 *    don't match.
 * 
 * Time complexity: O(nm), the size of the matrices.
 */

igraph_bool_t 
FUNCTION(igraph_matrix,all_le)(const TYPE(igraph_matrix) *lhs, 
			       const TYPE(igraph_matrix) *rhs) {
  return lhs->ncol==rhs->ncol && lhs->nrow==rhs->nrow &&
    FUNCTION(igraph_vector,all_le)(&lhs->data, &rhs->data);
}

/**
 * \function igraph_matrix_all_ge
 * \brief Are all elements greater or equal?
 * 
 * \param lhs The first matrix.
 * \param rhs The second matrix.
 * \return Positive integer (=true) if the elements in the \p lhs are all 
 *    greater than or equal to the corresponding elements in \p
 *    rhs. Returns \c 0 (=false) if the dimensions of the matrices
 *    don't match.
 * 
 * Time complexity: O(nm), the size of the matrices.
 */

igraph_bool_t 
FUNCTION(igraph_matrix,all_ge)(const TYPE(igraph_matrix) *lhs, 
			       const TYPE(igraph_matrix) *rhs) {
  return lhs->ncol==rhs->ncol && lhs->nrow==rhs->nrow &&
    FUNCTION(igraph_vector,all_ge)(&lhs->data, &rhs->data);
}

#endif

#ifndef NOTORDERED

/**
 * \function igraph_matrix_maxdifference
 * \brief Maximum absolute difference between two matrices.
 * 
 * Calculate the maximum absolute difference of two matrices. Both matrices
 * must be non-empty. If their dimensions differ then a warning is given and
 * the comparison is performed by vectors columnwise from both matrices.
 * The remaining elements in the larger vector are ignored.
 * \param m1 The first matrix.
 * \param m2 The second matrix.
 * \return The element with the largest absolute value in \c m1 - \c m2.
 * 
 * Time complexity: O(mn), the elements in the smaller matrix.
 */

BASE FUNCTION(igraph_matrix,maxdifference)(const TYPE(igraph_matrix) *m1,
					   const TYPE(igraph_matrix) *m2) {
  long int col1=FUNCTION(igraph_matrix,ncol)(m1);
  long int col2=FUNCTION(igraph_matrix,ncol)(m2);
  long int row1=FUNCTION(igraph_matrix,nrow)(m1);
  long int row2=FUNCTION(igraph_matrix,nrow)(m2);
  if (col1 != col2 || row1 != row2) {
    IGRAPH_WARNING("Comparing non-conformant matrices");
  }
  return FUNCTION(igraph_vector,maxdifference)(&m1->data, &m2->data);
}

#endif

/**
 * \function igraph_matrix_transpose
 * \brief Transpose a matrix.
 * 
 * Calculate the transpose of a matrix. Note that the function
 * reallocates the memory used for the matrix.
 * \param m The input (and output) matrix.
 * \return Error code.
 * 
 * Time complexity: O(mn), the number of elements in the matrix.
 */

int FUNCTION(igraph_matrix,transpose)(TYPE(igraph_matrix) *m) {
  long int nrow=m->nrow;
  long int ncol=m->ncol;
  if (nrow>1 && ncol>1) {
    TYPE(igraph_vector) newdata;
    long int i, size=nrow*ncol, mod=size-1;
    FUNCTION(igraph_vector,init)(&newdata, size);
    IGRAPH_FINALLY(FUNCTION(igraph_vector,destroy), &newdata);
    for (i=0; i<size; i++) {
      VECTOR(newdata)[i] = VECTOR(m->data)[ (i*nrow) % mod ];
    }
    VECTOR(newdata)[size-1]=VECTOR(m->data)[size-1];
    FUNCTION(igraph_vector,destroy)(&m->data);
    IGRAPH_FINALLY_CLEAN(1);
    m->data=newdata;
  }
  m->nrow=ncol;
  m->ncol=nrow;
  
  return 0;
}

/**
 * \function igraph_matrix_e
 * Extract an element from a matrix.
 * 
 * Use this if you need a function for some reason and cannot use the
 * \ref MATRIX macro. Note that no range checking is performed.
 * \param m The input matrix.
 * \param row The row index.
 * \param col The column index.
 * \return The element in the given row and column.
 * 
 * Time complexity: O(1).
 */

BASE FUNCTION(igraph_matrix,e)(const TYPE(igraph_matrix) *m, 
			       long int row, long int col) {
  return MATRIX(*m, row, col);
}

/**
 * \function igraph_matrix_e_ptr
 * Pointer to an element of a matrix.
 * 
 * The function returns a pointer to an element. No range checking is
 * performed. 
 * \param m The input matrix.
 * \param row The row index.
 * \param col The column index.
 * \return Pointer to the element in the given row and column.
 * 
 * Time complexity: O(1).
 */

BASE* FUNCTION(igraph_matrix,e_ptr)(const TYPE(igraph_matrix) *m,
				    long int row, long int col) {
  return &MATRIX(*m, row, col);
}

/**
 * \function igraph_matrix_set
 * Set an element.
 * 
 * Set an element of a matrix. No range checking is performed.
 * \param m The input matrix.
 * \param row The row index.
 * \param col The column index.
 * \param value The new value of the element.
 * 
 * Time complexity: O(1).
 */

void FUNCTION(igraph_matrix,set)(TYPE(igraph_matrix)* m, long int row, long int col,
				 BASE value) {
  MATRIX(*m, row, col) = value;
}

/**
 * \function igraph_matrix_fill
 * Fill with an element.
 * 
 * Set the matrix to a constant matrix.
 * \param m The input matrix.
 * \param e The element to set.
 * 
 * Time complexity: O(mn), the number of elements.
 */

void FUNCTION(igraph_matrix,fill)(TYPE(igraph_matrix) *m, BASE e) {
  FUNCTION(igraph_vector,fill)(&m->data, e);
}

/**
 * \function igraph_matrix_update
 * Update from another matrix.
 * 
 * This function replicates \p from in the matrix \p to. 
 * Note that \p to must be already initialized.
 * \param to The result matrix.
 * \param from The matrix to replicate; it is left unchanged.
 * \return Error code.
 * 
 * Time complexity: O(mn), the number of elements.
 */

int FUNCTION(igraph_matrix,update)(TYPE(igraph_matrix) *to, 
				   const TYPE(igraph_matrix) *from) {

  IGRAPH_CHECK(FUNCTION(igraph_matrix,resize)(to, from->nrow, from->ncol));
  FUNCTION(igraph_vector,update)(&to->data, &from->data);
  return 0;
}

/**
 * \function igraph_matrix_rbind
 * Combine two matrices rowwise.
 *
 * This function places the rows of \p from below the rows of \c to
 * and stores the result in \p to. The number of columns in the two
 * matrices must match.
 * \param to The upper matrix; the result is also stored here.
 * \param from The lower matrix. It is left unchanged.
 * \return Error code.
 * 
 * Time complexity: O(mn), the number of elements in the newly created
 * matrix.
 */

int FUNCTION(igraph_matrix,rbind)(TYPE(igraph_matrix) *to,
				  const TYPE(igraph_matrix) *from) {
  long int tocols=to->ncol, fromcols=from->ncol;
  long int torows=to->nrow, fromrows=from->nrow;
  long int offset, c, r, index, offset2;
  if (tocols != fromcols) {
    IGRAPH_ERROR("Cannot do rbind, number of columns do not match", IGRAPH_EINVAL);
  }
  
  IGRAPH_CHECK(FUNCTION(igraph_vector,resize)(&to->data, 
					      tocols * (fromrows+torows)));
  to->nrow += fromrows;

  offset=(tocols-1) * fromrows;
  index=tocols*torows-1;
  for (c=tocols-1; c>0; c--) {
    for (r=0; r<torows; r++, index--) {
      VECTOR(to->data)[index+offset] = VECTOR(to->data)[index];
    }
    offset -= fromrows;
  }
  
  offset=torows; offset2=0;
  for (c=0; c<tocols; c++) {
    memcpy(VECTOR(to->data)+offset, VECTOR(from->data)+offset2,
	   sizeof(BASE) * (size_t) fromrows);
    offset+=fromrows+torows;
    offset2+=fromrows;
  }
  return 0;
}

/**
 * \function igraph_matrix_cbind
 * Combine matrices columnwise.
 * 
 * This function places the columns of \p from on the right of \p to,
 * and stores the result in \p to.
 * \param to The left matrix; the result is stored here too.
 * \param from The right matrix. It is left unchanged.
 * \return Error code.
 * 
 * Time complexity: O(mn), the number of elements on the new matrix.
 */

int FUNCTION(igraph_matrix,cbind)(TYPE(igraph_matrix) *to,
				  const TYPE(igraph_matrix) *from) {
  
  long int tocols=to->ncol, fromcols=from->ncol;
  long int torows=to->nrow, fromrows=from->nrow;
  if (torows != fromrows) {
    IGRAPH_ERROR("Cannot do rbind, number of rows do not match", IGRAPH_EINVAL);
  }
  IGRAPH_CHECK(FUNCTION(igraph_matrix,resize)(to, torows, tocols+fromcols));
  FUNCTION(igraph_vector,copy_to)(&from->data, VECTOR(to->data)+tocols*torows);
  return 0;
}

/**
 * \function igraph_matrix_swap
 * Swap two matrices.
 * 
 * The contents of the two matrices will be swapped. They must have the
 * same dimensions.
 * \param m1 The first matrix.
 * \param m2 The second matrix.
 * \return Error code.
 * 
 * Time complexity: O(mn), the number of elements in the matrices.
 */

int FUNCTION(igraph_matrix,swap)(TYPE(igraph_matrix) *m1, TYPE(igraph_matrix) *m2) {
  if (m1->nrow != m2->nrow || m1->ncol != m2->ncol) {
    IGRAPH_ERROR("Cannot swap non-conformant matrices", IGRAPH_EINVAL);
  }
  return FUNCTION(igraph_vector,swap)(&m1->data, &m2->data);
}

/**
 * \function igraph_matrix_get_row
 * Extract a row.
 * 
 * Extract a row from a matrix and return it as a vector.
 * \param m The input matrix.
 * \param res Pointer to an initialized vector; it will be resized if
 *   needed.
 * \param index The index of the row to select.
 * \return Error code.
 * 
 * Time complexity: O(n), the number of columns in the matrix.
 */

int FUNCTION(igraph_matrix,get_row)(const TYPE(igraph_matrix) *m, 
				    TYPE(igraph_vector) *res, long int index) {
  long int rows=m->nrow, cols=m->ncol;
  long int i, j;
  
  if (index >= rows) {
    IGRAPH_ERROR("Index out of range for selecting matrix row", IGRAPH_EINVAL);
  }
  IGRAPH_CHECK(FUNCTION(igraph_vector,resize)(res, cols));
  
  for (i=index, j=0; j<cols; i+=rows, j++) {
    VECTOR(*res)[j] = VECTOR(m->data)[i];
  }
  return 0;
}

/**
 * \function igraph_matrix_set_row
 * Set a row from a vector.
 * 
 * Sets the elements of a row with the given vector. This has the effect of
 * setting row \c index to have the elements in the vector \c v. The length of
 * the vector and the number of columns in the matrix must match,
 * otherwise an error is triggered.
 * \param m The input matrix.
 * \param v The vector containing the new elements of the row.
 * \param index Index of the row to set.
 * \return Error code.
 * 
 * Time complexity: O(n), the number of columns in the matrix.
 */

int FUNCTION(igraph_matrix,set_row)(TYPE(igraph_matrix) *m,
				    const TYPE(igraph_vector) *v, long int index) {
  long int rows=m->nrow, cols=m->ncol;
  long int i, j;
  
  if (index >= rows) {
    IGRAPH_ERROR("Index out of range for selecting matrix row", IGRAPH_EINVAL);
  }
  if (FUNCTION(igraph_vector,size)(v) != cols) {
    IGRAPH_ERROR("Cannot set matrix row, invalid vector length", IGRAPH_EINVAL);
  }
  for (i=index, j=0; j<cols; i+=rows, j++) {
    VECTOR(m->data)[i]=VECTOR(*v)[j];
  }
  return 0;
}

/**
 * \function igraph_matrix_set_col
 * Set a column from a vector.
 *
 * Sets the elements of a column with the given vector. In effect, column
 * \c index will be set with elements from the vector \c v. The length of
 * the vector and the number of rows in the matrix must match,
 * otherwise an error is triggered.
 * \param m The input matrix.
 * \param v The vector containing the new elements of the column.
 * \param index Index of the column to set.
 * \return Error code.
 * 
 * Time complexity: O(m), the number of rows in the matrix.
 */

int FUNCTION(igraph_matrix,set_col)(TYPE(igraph_matrix) *m,
				    const TYPE(igraph_vector) *v, long int index) {
  long int rows=m->nrow, cols=m->ncol;
  long int i, j;
  
  if (index >= cols) {
    IGRAPH_ERROR("Index out of range for setting matrix column", IGRAPH_EINVAL);
  }
  if (FUNCTION(igraph_vector,size)(v) != rows) {
    IGRAPH_ERROR("Cannot set matrix column, invalid vector length", IGRAPH_EINVAL);
  }
  for (i=index*rows, j=0; j<rows; i++, j++) {
    VECTOR(m->data)[i]=VECTOR(*v)[j];
  }
  return 0;
}

/**
 * \function igraph_matrix_swap_rows
 * Swap two rows.
 * 
 * Swap two rows in the matrix.
 * \param m The input matrix.
 * \param i The index of the first row.
 * \param j The index of the second row.
 * \return Error code.
 *
 * Time complexity: O(n), the number of columns.
 */

int FUNCTION(igraph_matrix,swap_rows)(TYPE(igraph_matrix) *m, 
				      long int i, long int j) {
  long int ncol=m->ncol, nrow=m->nrow;
  long int n=nrow*ncol;
  long int index1, index2;
  if (i>=nrow || j>=nrow) {
    IGRAPH_ERROR("Cannot swap rows, index out of range", IGRAPH_EINVAL);
  }
  if (i==j) { return 0; }
  for (index1=i, index2=j; index1<n; index1+=nrow, index2+=nrow) {
    BASE tmp;
    tmp=VECTOR(m->data)[index1];
    VECTOR(m->data)[index1]=VECTOR(m->data)[index2];
    VECTOR(m->data)[index2]=tmp;
  }
  return 0;
}

/**
 * \function igraph_matrix_swap_cols
 * Swap two columns.
 * 
 * Swap two columns in the matrix.
 * \param m The input matrix.
 * \param i The index of the first column.
 * \param j The index of the second column.
 * \return Error code.
 *
 * Time complexity: O(m), the number of rows.
 */

int FUNCTION(igraph_matrix,swap_cols)(TYPE(igraph_matrix) *m, 
				      long int i, long int j) {
  long int ncol=m->ncol, nrow=m->nrow;
  long int k, index1, index2;
  if (i>=ncol || j >= ncol) {
    IGRAPH_ERROR("Cannot swap columns, index out of range", IGRAPH_EINVAL);
  }
  if (i==j) { return 0; }
  for (index1=i*nrow, index2=j*nrow, k=0; k<nrow; k++, index1++, index2++) {
    BASE tmp=VECTOR(m->data)[index1];
    VECTOR(m->data)[index1]=VECTOR(m->data)[index2];
    VECTOR(m->data)[index2]=tmp;
  }
  return 0;
}

/**
 * \function igraph_matrix_add_constant
 * Add a constant to every element.
 * 
 * \param m The input matrix.
 * \param plud The constant to add.
 * 
 * Time complexity: O(mn), the number of elements.
 */

void FUNCTION(igraph_matrix,add_constant)(TYPE(igraph_matrix) *m, BASE plus) {
  FUNCTION(igraph_vector,add_constant)(&m->data, plus);
}

/**
 * \function igraph_matrix_add
 * Add two matrices.
 * 
 * Add \p m2 to \p m1, and store the result in \p m1. The dimensions of the
 * matrices must match.
 * \param m1 The first matrix; the result will be stored here.
 * \param m2 The second matrix; it is left unchanged.
 * \return Error code.
 * 
 * Time complexity: O(mn), the number of elements.
 */

int FUNCTION(igraph_matrix,add)(TYPE(igraph_matrix) *m1, 
				const TYPE(igraph_matrix) *m2) {
  if (m1->nrow != m2->nrow || m1->ncol != m2->ncol) {
    IGRAPH_ERROR("Cannot add non-conformant matrices", IGRAPH_EINVAL);
  }
  return FUNCTION(igraph_vector,add)(&m1->data, &m2->data);
}

/**
 * \function igraph_matrix_sub
 * Difference of two matrices.
 * 
 * Subtract \p m2 from \p m1 and store the result in \p m1.
 * The dimensions of the two matrices must match.
 * \param m1 The first matrix; the result is stored here.
 * \param m2 The second matrix; it is left unchanged.
 * \return Error code.
 * 
 * Time complexity: O(mn), the number of elements.
 */

int FUNCTION(igraph_matrix,sub)(TYPE(igraph_matrix) *m1, 
				const TYPE(igraph_matrix) *m2) {
  if (m1->nrow != m2->nrow || m1->ncol != m2->ncol) {
    IGRAPH_ERROR("Cannot subtract non-conformant matrices", IGRAPH_EINVAL);
  }
  return FUNCTION(igraph_vector,sub)(&m1->data, &m2->data);
}

/**
 * \function igraph_matrix_mul_elements
 * Elementwise multiplication.
 * 
 * Multiply \p m1 by \p m2 elementwise and store the result in \p m1.
 * The dimensions of the two matrices must match.
 * \param m1 The first matrix; the result is stored here.
 * \param m2 The second matrix; it is left unchanged.
 * \return Error code.
 * 
 * Time complexity: O(mn), the number of elements.
 */

int FUNCTION(igraph_matrix,mul_elements)(TYPE(igraph_matrix) *m1, 
					 const TYPE(igraph_matrix) *m2) {
  if (m1->nrow != m2->nrow || m1->ncol != m2->ncol) {
    IGRAPH_ERROR("Cannot multiply non-conformant matrices", IGRAPH_EINVAL);
  }
  return FUNCTION(igraph_vector,mul)(&m1->data, &m2->data);
}

/**
 * \function igraph_matrix_div_elements
 * Elementwise division.
 * 
 * Divide \p m1 by \p m2 elementwise and store the result in \p m1.
 * The dimensions of the two matrices must match.
 * \param m1 The dividend. The result is store here.
 * \param m2 The divisor. It is left unchanged.
 * \return Error code.
 * 
 * Time complexity: O(mn), the number of elements.
 */

int FUNCTION(igraph_matrix,div_elements)(TYPE(igraph_matrix) *m1, 
					 const TYPE(igraph_matrix) *m2) {
  if (m1->nrow != m2->nrow || m1->ncol != m2->ncol) {
    IGRAPH_ERROR("Cannot divide non-conformant matrices", IGRAPH_EINVAL);
  }
  return FUNCTION(igraph_vector,div)(&m1->data, &m2->data);
}

#ifndef NOTORDERED

/**
 * \function igraph_matrix_min
 * Minimum element.
 * 
 * Returns the smallest element of a non-empty matrix. 
 * \param m The input matrix.
 * \return The smallest element.
 * 
 * Time complexity: O(mn), the number of elements.
 */

igraph_real_t FUNCTION(igraph_matrix,min)(const TYPE(igraph_matrix) *m) {
  return FUNCTION(igraph_vector,min)(&m->data);
}

/**
 * \function igraph_matrix_which_min
 * Indices of the minimum.
 * 
 * Gives the indices of the (first) smallest element in a non-empty
 * matrix.
 * \param m The matrix.
 * \param i Pointer to a <type>long int</type>. The row index of the 
 *   minimum is stored here.
 * \param j Pointer to a <type>long int</type>. The column index of
 *   the minimum is stored here.
 * \return Error code.
 * 
 * Time complexity: O(mn), the number of elements.
 */

int FUNCTION(igraph_matrix,which_min)(const TYPE(igraph_matrix) *m,
				      long int *i, long int *j) {
  long int vmin=FUNCTION(igraph_vector,which_min)(&m->data);
  *i = vmin % m->nrow;
  *j = vmin / m->nrow;
  return 0;
}

/**
 * \function igraph_matrix_which_max
 * Indices of the maximum.
 * 
 * Gives the indices of the (first) largest element in a non-empty
 * matrix.
 * \param m The matrix.
 * \param i Pointer to a <type>long int</type>. The row index of the 
 *   maximum is stored here.
 * \param j Pointer to a <type>long int</type>. The column index of
 *   the maximum is stored here.
 * \return Error code.
 * 
 * Time complexity: O(mn), the number of elements.
 */

int FUNCTION(igraph_matrix,which_max)(const TYPE(igraph_matrix) *m,
				      long int *i, long int *j) {
  long int vmax=FUNCTION(igraph_vector,which_max)(&m->data);
  *i = vmax % m->nrow;
  *j = vmax / m->nrow;
  return 0;
}

/**
 * \function igraph_matrix_minmax
 * Minimum and maximum
 * 
 * The maximum and minimum elements of a non-empty matrix.
 * \param m The input matrix.
 * \param min Pointer to a base type. The minimum is stored here.
 * \param max Pointer to a base type. The maximum is stored here.
 * \return Error code.
 * 
 * Time complexity: O(mn), the number of elements.
 */

int FUNCTION(igraph_matrix,minmax)(const TYPE(igraph_matrix) *m,
				   BASE *min, BASE *max) {
  return FUNCTION(igraph_vector,minmax)(&m->data, min, max);
}

/**
 * \function igraph_matrix_which_minmax
 * Indices of the minimum and maximum
 * 
 * Find the positions of the smallest and largest elements of a
 * non-empty matrix.
 * \param m The input matrix.
 * \param imin Pointer to a <type>long int</type>, the row index of
 *   the minimum is stored here.
 * \param jmin Pointer to a <type>long int</type>, the column index of
 *   the minimum is stored here.
 * \param imax Pointer to a <type>long int</type>, the row index of
 *   the maximum is stored here.
 * \param jmax Pointer to a <type>long int</type>, the column index of
 *   the maximum is stored here.
 * \return Error code.
 * 
 * Time complexity: O(mn), the number of elements.
 */

int FUNCTION(igraph_matrix,which_minmax)(const TYPE(igraph_matrix) *m,
					 long int *imin, long int *jmin,
					 long int *imax, long int *jmax) {
  long int vmin, vmax;
  FUNCTION(igraph_vector,which_minmax)(&m->data, &vmin, &vmax);
  *imin = vmin % m->nrow;
  *jmin = vmin / m->nrow;
  *imax = vmax % m->nrow;
  *jmax = vmax / m->nrow;
  return 0;
}

#endif

/**
 * \function igraph_matrix_isnull
 * Check for a null matrix.
 * 
 * Checks whether all elements are zero.
 * \param m The input matrix.
 * \return Boolean, \c TRUE is \p m contains only zeros and \c FALSE
 *   otherwise.
 *
 * Time complexity: O(mn), the number of elements.
 */

igraph_bool_t FUNCTION(igraph_matrix,isnull)(const TYPE(igraph_matrix) *m) {
  return FUNCTION(igraph_vector,isnull)(&m->data);
}

/**
 * \function igraph_matrix_empty
 * Check for an empty matrix.
 * 
 * It is possible to have a matrix with zero rows or zero columns, or
 * even both. This functions checks for these.
 * \param m The input matrix.
 * \return Boolean, \c TRUE if the matrix contains zero elements, and
 *    \c FALSE otherwise.
 * 
 * Time complexity: O(1).
 */

igraph_bool_t FUNCTION(igraph_matrix,empty)(const TYPE(igraph_matrix) *m) {
  return FUNCTION(igraph_vector,empty)(&m->data);
}

/**
 * \function igraph_matrix_is_symmetric
 * Check for symmetric matrix.
 * 
 * A non-square matrix is not symmetric by definition.
 * \param m The input matrix.
 * \return Boolean, \c TRUE if the matrix is square and symmetric, \c
 *    FALSE otherwise.
 * 
 * Time complexity: O(mn), the number of elements. O(1) for non-square
 * matrices. 
 */

igraph_bool_t FUNCTION(igraph_matrix,is_symmetric)(const TYPE(igraph_matrix) *m) {
  
  long int n=m->nrow;
  long int r,c;
  if (m->ncol != n) {
    return 0;
  }
  for (r=1; r<n; r++) {
    for (c=0; c<r; c++) {
      BASE a1=MATRIX(*m, r, c);
      BASE a2=MATRIX(*m, c, r);
#ifdef EQ      
      if (!EQ(a1,a2)) { return 0; }
#else
      if (a1 != a2) { return 0; }
#endif
    }
  }
  return 1;
}

/**
 * \function igraph_matrix_prod
 * Product of the elements.
 * 
 * Note this function can result in overflow easily, even for not too
 * big matrices.
 * \param m The input matrix.
 * \return The product of the elements.
 * 
 * Time complexity: O(mn), the number of elements.
 */

BASE FUNCTION(igraph_matrix,prod)(const TYPE(igraph_matrix) *m) {
  return FUNCTION(igraph_vector,prod)(&m->data);
}

/**
 * \function igraph_matrix_rowsum
 * Rowwise sum.
 * 
 * Calculate the sum of the elements in each row.
 * \param m The input matrix.
 * \param res Pointer to an initialized vector; the result is stored
 *   here. It will be resized if necessary.
 * \return Error code.
 * 
 * Time complexity: O(mn), the number of elements in the matrix.
 */

int FUNCTION(igraph_matrix,rowsum)(const TYPE(igraph_matrix) *m,
				   TYPE(igraph_vector) *res) {
  long int nrow=m->nrow, ncol=m->ncol;
  long int r, c;
  BASE sum;
  IGRAPH_CHECK(FUNCTION(igraph_vector,resize)(res, nrow));
  for (r=0; r<nrow; r++) {
    sum=ZERO;
    for (c=0; c<ncol; c++) {
#ifdef SUM
      SUM(sum, sum, MATRIX(*m, r, c));
#else
      sum += MATRIX(*m, r, c);
#endif
    }
    VECTOR(*res)[r]=sum;
  }
  return 0;
}

/**
 * \function igraph_matrix_colsum
 * Columnwise sum.
 * 
 * Calculate the sum of the elements in each column.
 * \param m The input matrix.
 * \param res Pointer to an initialized vector; the result is stored
 *   here. It will be resized if necessary.
 * \return Error code.
 * 
 * Time complexity: O(mn), the number of elements in the matrix.
 */

int FUNCTION(igraph_matrix,colsum)(const TYPE(igraph_matrix) *m,
				   TYPE(igraph_vector) *res) {
  long int nrow=m->nrow, ncol=m->ncol;
  long int r, c;
  BASE sum;
  IGRAPH_CHECK(FUNCTION(igraph_vector,resize)(res, ncol));
  for (c=0; c<ncol; c++) {
    sum=ZERO;
    for (r=0; r<nrow; r++) {
#ifdef SUM
      SUM(sum, sum, MATRIX(*m, r, c));
#else
      sum += MATRIX(*m, r, c);
#endif
    }
    VECTOR(*res)[c]=sum;
  }
  return 0;
}

/**
 * \function igraph_matrix_contains
 * Search for an element.
 * 
 * Search for the given element in the matrix.
 * \param m The input matrix.
 * \param e The element to search for.
 * \return Boolean, \c TRUE if the matrix contains \p e, \c FALSE
 * otherwise.
 * 
 * Time complexity: O(mn), the number of elements.
 */

igraph_bool_t FUNCTION(igraph_matrix,contains)(const TYPE(igraph_matrix) *m,
					       BASE e) {
  return FUNCTION(igraph_vector,contains)(&m->data, e);
}

/**
 * \function igraph_matrix_search
 * Search from a given position.
 * 
 * Search for an element in a matrix and start the search from the
 * given position. The search is performed columnwise.
 * \param m The input matrix.
 * \param from The position to search from, the positions are
 *    enumerated columnwise.
 * \param what The element to search for.
 * \param pos Pointer to a <type>long int</type>. If the element is
 *    found, then this is set to the position of its first appearance.
 * \param row Pointer to a <type>long int</type>. If the element is
 *    found, then this is set to its row index.
 * \param col Pointer to a <type>long int</type>. If the element is
 *    found, then this is set to its column index.
 * \return Boolean, \c TRUE if the element is found, \c FALSE
 *    otherwise.
 * 
 * Time complexity: O(mn), the number of elements.
 */

igraph_bool_t FUNCTION(igraph_matrix,search)(const TYPE(igraph_matrix) *m,
					     long int from, BASE what, 
					     long int *pos, 
					     long int *row, long int *col) {
  igraph_bool_t find=FUNCTION(igraph_vector,search)(&m->data, from, what, pos);
  if (find) {
    *row = *pos % m->nrow;
    *col = *pos / m->nrow;
  }
  return find;
}

/**
 * \function igraph_matrix_remove_row
 * Remove a row.
 * 
 * A row is removed from the matrix. 
 * \param m The input matrix.
 * \param row The index of the row to remove.
 * \return Error code.
 * 
 * Time complexity: O(mn), the number of elements in the matrix.
 */

int FUNCTION(igraph_matrix,remove_row)(TYPE(igraph_matrix) *m, long int row) {
  
  long int c, r, index=row+1, leap=1, n=m->nrow * m->ncol;
  if (row >= m->nrow) {
    IGRAPH_ERROR("Cannot remove row, index out of range", IGRAPH_EINVAL);
  }  
  
  for (c=0; c<m->ncol; c++) {    
    for (r=0; r<m->nrow-1 && index < n; r++) {
      VECTOR(m->data)[index-leap] = VECTOR(m->data)[index];
      index++;
    }
    leap++;
    index++;
  }
  m->nrow--;
  FUNCTION(igraph_vector,resize)(&m->data, m->nrow * m->ncol);
  return 0;
}

/**
 * \function igraph_matrix_select_cols
 * \brief Select some columns of a matrix.
 * 
 * This function selects some columns of a matrix and returns them in a
 * new matrix. The result matrix should be initialized before calling
 * the function.
 * \param m The input matrix.
 * \param res The result matrix. It should be initialized and will be
 *    resized as needed.
 * \param cols Vector; it contains the column indices (starting with
 *    zero) to extract. Note that no range checking is performed.
 * \return Error code.
 * 
 * Time complexity: O(nm), n is the number of rows, m the number of
 * columns of the result matrix.
 */

int FUNCTION(igraph_matrix,select_cols)(const TYPE(igraph_matrix) *m,
					TYPE(igraph_matrix) *res, 
					const igraph_vector_t *cols) {
  long int ncols=igraph_vector_size(cols);
  long int nrows=m->nrow;
  long int i, j;
  
  IGRAPH_CHECK(FUNCTION(igraph_matrix,resize)(res, nrows, ncols));
  for (i=0; i<nrows; i++) {
    for (j=0; j<ncols; j++) {
      MATRIX(*res, i, j) = MATRIX(*m, i, (long int)VECTOR(*cols)[j]);
    }
  }
  return 0;
}

#ifdef OUT_FORMAT

#ifndef USING_R
int FUNCTION(igraph_matrix,print)(const TYPE(igraph_matrix) *m) {

  long int nr=FUNCTION(igraph_matrix,nrow)(m);
  long int nc=FUNCTION(igraph_matrix,ncol)(m);
  long int i, j;
  for (i=0; i<nr; i++) {
    for (j=0; j<nc; j++) {
      if (j!=0) { putchar(' '); }
      printf(OUT_FORMAT, MATRIX(*m, i, j));
    }
    printf("\n");
  }
  
  return 0;
}

int FUNCTION(igraph_matrix,printf)(const TYPE(igraph_matrix) *m,
																	 const char *format) {
  long int nr=FUNCTION(igraph_matrix,nrow)(m);
  long int nc=FUNCTION(igraph_matrix,ncol)(m);
  long int i, j;
  for (i=0; i<nr; i++) {
    for (j=0; j<nc; j++) {
      if (j!=0) { putchar(' '); }
      printf(format, MATRIX(*m, i, j));
    }
    printf("\n");
  }

  return 0;
}

#endif

int FUNCTION(igraph_matrix,fprint)(const TYPE(igraph_matrix) *m,
				   FILE *file) {

  long int nr=FUNCTION(igraph_matrix,nrow)(m);
  long int nc=FUNCTION(igraph_matrix,ncol)(m);
  long int i, j;
  for (i=0; i<nr; i++) {
    for (j=0; j<nc; j++) {
      if (j!=0) { fputc(' ', file); }
      fprintf(file, OUT_FORMAT, MATRIX(*m, i, j));
    }
    fprintf(file, "\n");
  }
  
  return 0;
}

#endif
